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ABSTRACT 

We present illustrative models for the UV to millimeter emission of starburst galaxies 
which are treated as an ensemble of optically thick giant molecular clouds (GMCs) cen- 
trally illuminated by recently formed stars. The models follow the evolution of the GMCs 
due to the ionization-induced expansion of the HII regions and the evolution of the stellar 
population within the GMC according to the Bruzual & Chariot stellar population syn- 
thesis models. The effect of transiently heated dust grains/PAHs to the radiative transfer, 
as well as multiple scattering, is taken into account. 

The expansion of the HII regions and the formation of a narrow neutral shell naturally 
explains why the emission from PAHs dominates over that from hot dust in the near to 
mid-IR, an emerging characteristic of the infrared spectra of starburst galaxies. 

The models allow us to relate the observed properties of a galaxy to its age and star 
formation history. We find that exponentially decaying 10 7 — 10 8 yrs old bursts can explain 
the IRAS colours of starburst galaxies. The models are also shown to account satisfactorily 
for the multiwavelength data on the prototypical starburst galaxy M82 and NGC6090, a 
starburst galaxy recently observed by ISO. In M82 we find evidence for two bursts separated 
by 10 7 yrs. In NGC6090 we find that at least part of the far-IR excess may be due to the age 
of the burst ( 6.4 x 10 7 yrs). We also make predictions about the evolution of the luminosity 
of starbursts at different wavelengths which indicate that far-IR surveys may preferentially 
detect older starbursts than mid-IR surveys. 
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1 INTRODUCTION 

It now looks increasingly clear that the long-sought era of 
galaxy formation is becoming accessible with a number of ob- 
servational methods for close scrutiny. These range from op- 
tical/ultraviolet studies (Steidel & Hamilton 1992, Lilly et al 
1996, Madau et al 1996), mid-infrared (Rowan- Robinson et 
al 1997 and references therein) and sub-millimeter surveys 
(Hughes et al 1998). A number of authors have discussed the 
history of star formation in the Universe which appears to have 
peaked at a redshift of about 1-3. Estimates like this are, of 
course, model-dependent as they involve the converion from an 
observed (usually monochromatic) luminosity density to a star- 
formation rate. In particular, they depend on the extent and 
geometry of dust obscuration especially at the shorter wave- 
lengths. 

As is well known a large fraction of the power emitted 
by galaxies (ranging from about 30% in normal galaxies to 
almost 100% in actively star forming galaxies or starbursts) 
lies in the infrared part of the spectrum as a result of repro- 
cessing of starlight by dust. Extensive infrared observations 
of galaxies are therefore necessary in order to describe fully 
their energy outnut. With the advent of the Infrared Soace 



local Universe have been observed with unprecedented detail 
in the infrared. They invariably display a variety of absorp- 
tion/emission features due to dust/molecules. It is clear that 
radiative transfer calculations in dusty media will be useful for 
the interpretation of these observations and a development of 
a better understanding of the origin of the infrared luminosity 
of galaxies. 

Radiative transfer models for the infrared emission of star- 
burst galaxies have been presented before by Rowan-Robinson 
& Crawford (1989; hereafter RRC), Rowan- Robinson & Efs- 
tathiou (1993; hereafter RRE), Kriigel & Siebenmorgen (1994; 
hereafter KS). These models used state of the art codes for 
calculating the transfer of radiation in dusty media and in- 
corporated a model for the composition and size distribution 
of grains in the interstellar medium. KS additionally included 
the effect of transiently heated grains/PAHs and considered 
the local change of dust temperatures in hot spots around 
luminous OB stars. The basic assumption of previous star- 
burst models is that a starburst is made up of an ensemble 
of compact HII regions similar to those found in our galaxy. 
A number of authors (e.g. Rowan-Robinson 1980, Kriigel & 
Walmslev 1984, Churchwell, Wolfire & Wood 1990, Efstathiou 
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different samples of such HII regions and concluded that they 
could be modelled adequately by spherically symmetric dust 
clouds surrounding young massive stars. 

All of the above models assume that the cloud ensemble in 
the starburst consists of a number of identical systems. In this 
paper we present illustrative models for the evolution of giant 
molecular clouds (GMCs) induced by massive star formation 
at their centers and calculate their infrared spectra. In the 
radiative transfer code we use, the effect of transiently heated 
particles/PAHs as well as classical grains is included. We also 
follow the evolution of the stellar populations with the models 
of Bruzual & Chariot (1995). This approach allows us to relate 
the observed properties of a starburst to its age and its star 
formation history. We illustrate these models by comparison 
with the IRAS colours of starburst galaxies as a class and with 
the multiwavelength data on M82 and NGC6090. 



2 A NEW STARBURST MODEL 

The basic assumption of our model for a starburst galaxy is 
that star formation takes place primarily within optically thick 
molecular clouds. This is supported by an array of observa- 
tional studies (e.g. Elmegreen 1985) which show that molecular 
clouds are associated with young stars. 

The mass of molecular clouds in our Galaxy ranges be- 
tween 10 2 — 10 7 M Q with a mass distribution approximately 
following M" 1 - 5 * - 1 (Dame et al 1983, Solomon et al 1987). It 
follows that about 70% of molecular mass is associated with 
GMCs more massive than 10 6 M Q . If the mass distribution in a 
starburst follows a similar form as the Galactic one, and GMCs 
are at least as likely to be the sites of massive stars as the less 
massive molecular clouds, then we would expect the bulk of the 
luminosity of the starburst to arise from GMCs with a fairly 
narrow range of mass. In fact low-mass molecular clouds, like 
the nearby Taurus and Ophiuchus clouds, are known to form 
predominantly low-mass stars. 

An indirect argument that massive star formation in 
galaxies takes place in GMCs with roughly a power-law mass 
distribution comes from HST observations of starburst galaxies 
(O'Connell et al 1995, Meurer et al 1995) which reveal a pop- 
ulation of super star clusters following a luminosity function 
of a power-law form 4>(L) oc L~ 2 . A similar slope is found for 
systems of young clusters in other galaxies (e.g. Whitmore & 
Schweizer 1995). This luminosity function is quite unlike that 
of globular clusters although it has been suggested that it could 
evolve into one (Meurer et al). In our model these star clus- 
ters represent the evolved counterparts of star clusters forming 
within GMCs, with the mass spectrum given above, after they 
have dispersed their nascent molecular clouds. It is not clear 
what the ages of these clusters are because of the known age- 
reddening degeneracy but in the case of M82 their proximity 
to the 2/im nucleus (see section 4) points towards an age > 
lOMyrs. 

In previous studies we assumed that a starburst is an en- 
semble of identical star forming complexes, which we approxi- 
mate as spherical. Here we refine this model by considering the 
fact that, given the starburst takes place over a finite period of 
time, the star-forming complexes that constitute it are bound 
to be at different evolutionary stages. We therefore use a sim- 
ple evolutionary model for HII regions to construct a family of 
models that predict the infrared spectrum of HII regions as a 
function of their age. 



2.1 The dusty HII region phase 

Star formation takes place primarily in the dense cores of 
GMCs. The details of the physical processes involved are not 
well understood yet. The efficiency of star formation (or the gas 
consumption rate) ranges from about 1% in late type spirals 
to 60% or more in starburst galaxies (Kennicutt 1998). 

Once massive stars form they can ionize the surrounding 
medium and inhibit further star formation. This in itself, how- 
ever, cannot explain the low star formation efficiency in disk 
galaxies because it does not explain what determines the num- 
ber of massive young stars formed initially. 

The evolution of GMCs is determined by ionization in- 
duced expansion in the early stages (t < 10 7 yrs) and later 
by stellar winds and supernova explosions. The latter eventu- 
ally disperse the molecular clouds on timescales of a few times 
10 7 years. Current scenaria of self-propagating star formation 
(e.g. Larson 1988, Tenorio-Tagle & Bodenheimer 1988, here- 
after TTB) hold that molecular cloud formation by a number 
of mechanisms occurs on a timescale of 10 s yrs, to complete 
the cycle. 

The evolution of HII regions due to ionization has been the 
subject of extensive study both analytically and numerically 
for a broad variety of circumstances (spherically symmetric 
case, two-dimensional solutions leading to so-called champagne 
flows) and is considered to be generally well understood (for a 
review see Franco et al 1990, TTB). 

The presence of a number of massive stars in the centre 
of a CMC producing numerous ionizing photons leads to the 
formation, on an extremely short timescale, of the classical HII 
region with initial Stromgren radius Rs- The pressure in the 
ionized region, being orders of magnitudes higher than in the 
rest of the GMC, drives the expansion of the HII region into 
the surrounding medium. The details of the expansion depend 
on a number of factors including the density distribution in 
the medium. Assuming a constant density medium and spher- 
ical symmetry, the evolution of the HII region radius can be 
described by (Spitzer 1978) 



where a is the sound speed in the ionized gas. t is the time since 
the formation of the initial Stromgren sphere. However, since 
the timescale for the latter is extremely short we will assume 
that t is the time from the onset of star formation (assumed 
to take place instantaneously) in the GMC. 

The expansion can be broken down into a number of 
phases, marked by the time when massive stars first move 
off the main sequence and some recombination occurs (t ~ 
3 — 4 x 10 6 yrs), a phase of further ionization, and finally to- 
tal recombination (t w 10 7 yrs) as F, diminishes. However, 
throughout expansion and until the supernovae ejecta finally 
reach the expansion front caused by ionization, equation (1) 
can be considered to give a reasonable description of the evo- 
lution, a fact confirmed by numerical simulations (TTB). 

In the simple case of a core of constant density n c crn~ 3 
the initial Stromgren radius Rs generated by a compact star 
cluster producing F» ionizing photons is given by (Franco et 
al 1990), 

Rs = A - 9 ( sxlO^-O " ( 2xlL-3 ) 3 PC (2) 

The average densities of GMCs in our Galaxy are in the 
range 10 to 10 2 cm" 3 (Dame et al 1983) but their cores, where 
most of the stars form, have densities three or more orders of 
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actively star-forming galaxies in accordance with the Schmidt 
law (Kennicutt 1998). Given the lack of knowledge about n c , 
its density distribution, the effect of dust within the HI1 region 
and our assumption that all star formation in a given GMC 
occurs instantaneously on defining Rs, equation (2) provides 
only a rough estimate of Rs, probably an upper limit. 

We use the tables of Bruzual and Chariot (1995) to derive 
the number of Lyman continuum photons F*. For a Salpeter 
IMF and stars in the mass range 0.1 — 125M©, 
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where r\ is the star formation efficiency. 

The molecular cloud radius r 2 is related to its mass Mqmc 
and average density n av (initially assumed to be uniform) by 
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At time t = 0, the GMC is essentially divided into three 
zones. Zone A (inside the dust sublimation radius n) consists 
of ionized gas and the stellar cluster which is approximated as 
a point source. In zone B ( n < r < Rs) we have ionized gas 
and dust, whereas zone C ( Rs < r < n) is the dusty neutral 
zone. 

As the expansion gets under way, the shock wave at R(t), 
leading the ionization front at Ri(t), accumulates neutral gas 
between Ri and R. Inside Ri the mass density of the ions pi 
follows (Franco ct al ) 



Pi(t) = p no (^) 
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(5) 



where p„ is the initial mass density of the neutral material 
assumed to be equal to the initial mass density of the ions. 
The density inside the HII region is probably constant to a 
good approximation as the expansion is subsonic (Franco et al 
)• 

The density enhancement and the separation of Ri and 
R cannot be followed analytically. Instead we assume that the 
neutral gas accumulated in the shock is spread over the entire 
neutral cloud, and we use the principle of conservation of mass 
to calculate the (uniform) density of neutral material p n (t) 
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Assuming the standard conversion from gas column den- 
sity to extinction (Savage & Mathis 1979), the visual extinction 
to the center of the cloud at t = is given by, 



tv = 50 C 



( 
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(7) 



where £ is the metallicity with respect to solar. 

To estimate the evolution of the optical depth of the cloud 
we assume that the number density of grains scales in the same 
way as pi and p n to get 



T. = nU^(i--) +/ d ^(--^)] 

Pn r 2 p no r 2 r 2 



(8) 



where fd takes account of a possible depletion of dust in the 
HII region as a result of shocks etc. 

We will assume that the shock advances into the neutral 
part of the GMC until the swept up material is half of the 
original GMC mass. For our assumed uniform densitv medium 



2.2 The supernova phase 

The state of the GMC at ~ 10 7 years is likely to be a narrow 
neutral shell still expanding, because of the momentum it ac- 
quired in the HII phase, at about 10 km s _1 . Over the next 
3 — 4 x 10 7 years supernova explosions are likely to be the main 
factor that will influence the evolution of the GMCs eventually 
leading to their dispersal and return of most of the gas to the 
HI phase. 

This phase has also been extensively studied analytically 
and numerically, but unfortunately mainly for the 'standard' 
density case of 1 cm~ 3 (TTB and references therein). Most 
relevant for our purposes is the work of McCray and Kafatos 
(1987) who presented an analytical solution for the case of 
multi-supernova explosions in an OB association. Numerical 
solutions have shown that this solution is an adequate approx- 
imation if the supernova explosions are frequent as in our case. 
The sequence of discrete supernova explosions in the model is 
replaced by a scaled-up stellar wind solution. The radius of the 
swept up shell before the hot interior begins to cool at time t c 
is given by 
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where N* is the total number of stars in the association (with 
M > 7Mq) that are destined to explode as supernovae , E51 is 
the energy per supernova (assumed constant) in units of 10 51 
ergs, and n is the density of the medium, again assumed to be 
constant. 

After time t c , given by 

(10) 



4 x 10 6 r 1 " 5 (iV.£5i) - 3 n- -V 



the shell expands as 
Rsn = Rc(t/t c ) 1/4 



(11) 



where R c is the radius reached at t c . 

For a Salpeter IMF and stars in the mass range 0.1 — 
125M© the Bruzal & Chariot models predict 8.92 x 10~ 3 su- 
pernova explosions per 1 M of stellar mass. So, our adopted 
parameters for the GMCs in a starburst (Mgmc = 10 7 M Q , 
77 = 0.25; see section 3) imply iV* = 2.23 x 10 4 . This, in 
turn, implies that the supernova shell will remain inside the 
shell formed by ionization at t = 10 7 yrs if the density of the 
medium within which the supernovae are exploding is higher 
than 10 3 cm -3 . This is not unreasonable given the contribu- 
tions from stellar winds, failed 'cores' etc. in a compact cluster 
of this mass. The velocity of the supernova shell is also pre- 
dicted to drop well below the expansion velocity of the neutral 
shell by 10 7 yrs. 

In this paper we will therefore assume that the supernova 
shell remains 'trapped' inside the neutral shell which continues 
to expand at 10 km s _1 until t = 4x 10 7 yrs. We discuss whether 
this assumption is inconsistent with infrared observations of 
galaxies later in the paper. 

Our plan in this paper is to use the simple evolutionary 
scenario outlined in the last two sections and our radiative 
transfer code in dusty media to predict the infrared properties 
of evolving HII regions. A further ingredient of our models will 
be the evolutionary synthesis models of Bruzual & Chariot 
(1995) which will provide the spectral energy distribution of 
the evolving stellar population powering the HII region. 

The family of evolving HII region models will then form 
the basis of our starburst models under our assumption that 
the latter are an ensemble of giant HII regions at different 
evolutionarv stages. 
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our models we have in this paper confined our attention to 
the spherically symmetric case. Clearly, highly non-spherical 
geometries can arise especially at the latter stages of the evo- 
lution and we plan to explore these situations in future studies. 



2.3 Dust model 

The model we have assumed for the absorption/emission prop- 
erties of the dust is an extension of the 'classical' grain model 
(e.g. Mathis 1990 and references therein) to take into account 
the effects of small grains and molecules. The model is de- 
scribed in detail in Siebenmorgen & Kriigel (1992; hereafter 
SK92) where it has been shown to account satisfactorily for 
the emitted spectra of dust in a number of environments (so- 
lar neighbourhood, planetary nebulae, star-forming regions). 

The model assumes three populations of grains and aims 
to fit the average interstellar extinction curve subject to the 
constraints imposed by abundances of heavy elements that are 
in the solid state etc. The first population of grains consists of 
the large grains (assumed to have optical properties as the "as- 
tronomical" silicates of Draine & Lee 1984 and of amorphous 
carbon by Edoh 1983). These grains provide the bulk of the 
emission/absorption at long wavelengths and are responsible 
for the silicate resonances at 9.7 and 18pim. The large grains 
are assumed to have a power-law size distribution (n(a) cx a~ 9 , 
q = 3.5, 100A< a <2500A). 

The second population of grains consists of small graphite 
particles (n(a) oc a' ~ q , q = 4, 10A< a <100A). These grams 
are responsible for the 2175A bump (Draine 1989) and emit 
primarily in mid-IR wavelengths. Because of their small size 
they show dramatic fluctuations in temperature. Their emis- 
sion is calculated with the method of Siebenmorgen et al (1992) 
which is a faster (but still fairly accurate) alternative to the 
treatment of Guhathakurta & Draine (1989). 

The third grain population is composed of PAHs that are 
now widely believed to be the carriers of the infrared features 
at 3.3, 6.2, 7.7, 8.6 and 11.3/ttm (Puget & Leger 1989, Alla- 
mandola et al 1989). The SK92 model assumes two compo- 
nents for the PAHs: a single molecule component made up of 
25 carbon atoms and a cluster component made up of 10-20 
molecules. A ratio of carbon atoms in PAHs to H atoms in the 
ISM of 3 x 10~° is assumed (which is about 10% of the total 
carbon abundance thought to reside in grains/molecules). The 
PAH clusters dominate at longer (mid-IR) wavelengths and 
are largely responsible for the quasi-continuum underlying the 
features (Desert, Boulanger & Puget 1990). The PAHs absorb 
mainly in the Far UV and to a lesser extent in the visible 
(SK92). 

The PAHs are thought to be dehydrogenated depending 
on their environment, so a degree of dehydrogenation oih/c 
(defined as the ratio of hydrogen to carbon atoms) needs to 
be introduced, olh/c varies from about 0.4 for the solar neigh- 
bourhood to 0.1 or less for star forming regions (SK92). As the 
emission of some of the features (3.3, 8.6 and 11.3/im) is pro- 
portional to the number of hydrogen atoms and of others (6.2 
and 7.7/xm) proportional to the number of carbon atoms, a H / C 
can have an effect on the relative strengths of the features. 

In regions of high radiation intensity (such as those to be 
found in regions of massive star formation or AGN), the small 
graphites and PAHs are thought to be underabundant with 
respect to the ISM by factors of 10 or more (SK92, Rowan- 
Robinson 1992 and references therein). This will result in a 
flattening of the extinction curve in the UV and an elimina- 
tion of the 2175A feature. The latter has been shown to be 



(Gordon, Calzetti & Witt 1997). For the models presented in 
this paper we assume that all the grains are depleted by a fac- 
tor of 5 (i.e. f d — 0.2) inside R. We further assume that the 
clusters are made of 500 C atoms and contain 90% of the total 
C abundance in PAHs. Note that the PAHs abundance is fur- 
ther reduced because of the photo-destruction mechanism that 
is self-consistently applied in the code (Siebenmorgen 1993). 

2.4 Radiative transfer model 

The method of solution of the radiative transfer problem in 
dusty media containing transiently heated grains is described 
in Efstathiou & Siebenmorgen (1999). The method of obtain- 
ing the intensity distribution at any point in the cloud and 
hence iterating for the temperature of each of the large grains 
is that used by Efstathiou & Rowan-Robinson (1990, 1995; 
hereafter ER90, ER95 respectively). The emission of the tran- 
siently heated particles is calculated according to the method 
of Siebenmorgen et al. (1992). Proper treatment of the pho- 
todestruction of the PAHs and the sublimation of the large 
grains, at the inner part of the cloud, is taken into account. 



3 EVOLVING HII REGION MODELS 

There are basically three free parameters in our model (Mgmc, 
j), and n av ) which we fix by relating to observational con- 
straints where available. 

Kennicutt (1998) finds the median rate of gas consumption 
in starburst galaxies per 10 8 years to be ~ 30%. Colbert et 
al (1998) estimate from ISO spectroscopy that the mass in 
stars in the M82 starburst is 0.5 — 1.3 x 1O 8 M . Assuming the 
molecular gas mass of ~ 2 x 10 8 Mq estimated from a number 
of studies (e.g. Hughes et al 1994) is associated with the GMCs 
that formed those stars, the implied gas consumption rate is 
0.2-0.39. In this paper we assume r\ = 0.25. 

For Mgmc we choose 1O 7 M , a mass close to the upper 
limit of the galactic mass spectrum, as the emission is likely 
to be dominated by such GMCs (section 2). As a check, we 
can compare the luminosity expected from one of these GMCs 
{Mgmc = 1O 7 M ,?7 = 0.25) with that of the 10^m knots in 
M82. Assuming they have the same spectrum as the starburst 
as a whole, each of these knots has a luminosity of ~ 2 x 10 L 
(Telesco & Gezari 1992). By comparison, the models of Bruzual 
& Chariot (0.1 < M < 125M ) predict L ~ 3.5 x 1O 8 L for a 
10 7 years old instantaneous burst. 

Probably the most uncertain free parameter in our mod- 
els is the initial CMC average density n av and the core den- 
sity n c . We assume n av = 300cm -3 which is within the range 
of GMC densities found in the centre of the galaxy (Glisten 
1989). There is evidence (e.g. RRE, Downes & Solomon 1998) 
that higher gas densities are to be found in the more extreme 
starbursts powering ultraluminous infrared galaxies. The core 
density is assumed to be 2 x 10 3 cm -3 . 

In Figure 1 we plot the density, temperature and spec- 
tral energy distributions of a GMC with assumed parameters 
Mgmc = 1O 7 M ,?7 = 0.25, and n av = 300cm -3 . These pa- 
rameters are used for the remainder of this paper. The density 
at the boundary between the neutral and ionized region has 
been smoothed in order to be able to handle it numerically. 
The shading in the same diagram represents the spread in tem- 
perature between different grain species. Note that the spread 
is much smaller in the neutral region because it is more op- 
tically thick (see Efstathiou & Rowan-Robinson 1994, Kriigel 
& Walmslev 1984 for a discussion of such radiative transfer 
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Figure 1. Temperature and density distributions (left) and spectral energy distributions (right), of the GMCs at two representative ages of 
the HII region phase. The shaded regions show the spread in temperature attained by the different grain species. The density and vSv scale 
is arbitrary. 



by 10 years the ionization front has compressed the dust to a 
very narrow shell (n/r2 ~ 0.8). 

The SEDs of the GMCs vary significantly with age. In the 
early stages of the cloud's evolution its SED is predicted to be 
warmer and show little signs of PAH emission. This is partly 
because the stellar population is younger and the radiation 
field stronger. The main contributing factor though is that 
there is more hot dust inside the Stromgren sphere than at 
later times. The weakness of the PAH features is partly due 
to the stronger near to mid-IR continuum emission from large 
grains and to the higher degree of photodestruction because of 
the stronger and harsher radiation field. 

By lOMyrs the mid-IR spectrum is dominated by the 
PAHs features and quasi-continuum and shows the character- 
istic shape that is observed in the spectra of starburst galaxies 
(Willner et al 1977, Acosta-Pulido et al 1996). Note that while 
the general trend is for the dust shell to cool off with time, 
the expansion of the HII region introduces a subtle effect that 
somewhat counteracts this in the HII phase. As the neutral 
shell is pushed out and the density inside R declines, the to- 
tal optical depth to the stellar cluster (eqn 8) decreases. This 
means that the neutral shell actually heats up a little bit. This 
effect is better demonstrated by the evolution of the IRAS 
colours, plotted in Figure 2. 

At t > 20Myrs the peak of the SEDs of the shells shifts 
to longer wavelengths and become remarkably similar to those 
of the diffuse medium and cirrus clouds. 

The effect of changing the IMF from Salpeter to Miller- 
Scalo is to reduce F, by about 80% and Rs by about 20%. So 
the effect on the overall SEDs is small. 



4 EVOLVING STARBURST MODELS 

To synthesize the spectral energy distribution from a burst of 
star formation from those of individual GMCs, let us assume 
that at time t after the onset of the starburst the star forma- 
tion rate (or in our case the number of GMCs forming stars 
instantaneously with efficiency rf) is Ai»(t). If we further as- 
sume that the ensemble of clouds is optically thin, i.e. they 
don't shadow each other, then the flux from the burst is given 
by 

F v {t)= I M*(t')S„{t-t')dt' (12) 
J o 

where S v (t —t') is the flux from a GMC t — t' years after the 
onset of star formation at its centre. 

A useful parameterization for the star formation rate in 
a starburst which has been extensively used (e.g. Rieke et al 
1980, Genzel et al 1998) is that of exponential decay, M*(t') ~ 
e~*/ T , where r is some time constant. In the limit of large r 
this approximates a scenario of constant star formation history. 
Under this assumption eqn (12) reduces to 

F v (t) = M,(0)e~ t/T ( e* ' /t S v (t')dt' (13) 
Jo 

In general, some stellar light will escape the GMC without 
being absorbed by the dust shell, a situation we can't address 
exactly in our present spherically symmetric model. To account 
for this, but only in the models fitted to M82 and NGC6090, 
S v is corrected for this effect by allowing a fraction f(t) of 
the starlight to leak out. The emission from the dust is cor- 
respondingly reduced to maintain flux conservation. The leak- 
ing starlight may suffer further absorption with the radiation 
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Figure 2. IRAS colour-colour diagrams of starburst galaxies from the sample of Rowan-Robinson & Crawford (1989) with models superim- 
posed. The dotted line shows the colours of the sequence of GMCs. The ages indicated by open squares are (clockwise on the 25/12 versus 
60/25 colour plot) 0, 1.7, 6.6, 10, 17, 26, 32, 37, 45, 57 and 72Myrs. The solid line shows the colours of a r = 20Myrs starburst, calculated 
using equation (12), at the same ages. The symbol D denotes the position of the 'Disc' component of RRC. To obtain the model colours the 
model SEDs were convolved with the IRAS filter response functions as explained in the IRAS Explanatory Supplement. 



number of possibilities here. The starlight may be absorbed 
locally (i.e. the GMC has some kind of anisotropic density dis- 
tribution, or optically thin holes), or alternatively the starlight 
escapes the GMC and is absorbed by another GMC or the dif- 
fuse medium (KS). In the first case the best approximation 
would probably be to assume that the starlight is reprocessed 
to the infrared with the same SED as that of the GMC it- 
self. In the second case it would be better to assume that the 
starlight is reprocessed according to the SED of the more nu- 
merous GMCs with ages of about 40-100Myrs. In this paper, 
in particular in the models fitted to M82 and NGC6090 in 
section 6, we consider only the first case. More details on the 
approximation used in those fits can be found in Appendix A. 



5 COMPARISON WITH IRAS DATA 

The IRAS data on galaxies have for the last decade or so 
provided the yardstick by which theoretical models of their 
infrared emission should be measured. While more extensive 



other projects (e.g. SCUBA), which enable more detailed com- 
parisons (and we illustrate this in the next section for two well 
studied objects), it is still instructive to compare our models 
with the IRAS colours. 

RRC proposed that the IRAS colours of galaxies can be 
understood in terms of 3 components: a 'disc' component with 
infrared properties similar to those of cirrus clouds in our 
Galaxy, a starburst component with properties similar to those 
of compact galactic HII regions and an AGN component, which 
is associated with the dusty torus that is now an integral part 
of the standard AGN model. The latter component is not dis- 
cussed further in this paper but instead we concentrate on the 
colours of starburst and Disc galaxies. 

In Figure 2 we plot the IRAS colours of the starburst 
galaxies in the sample of RRC (selected to have good qual- 
ity fluxes in all four IRAS bands) and indicate the position of 
their disc component (D). Also indicated in Figure 2 are the 
positions of M82 and NGC6090 on the colour-colour diagrams. 
Disc galaxies cover the Dart of the colour-colour diagram be- 
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Figure 3. Spectral energy distributions of a r = 20Myrs starburst at different ages. The assumed GMC parameters are Mqmc = 10 7 Mq, r\ = 
0.25, riav = 300cra~ 3 . The vertical scale is arbitrary but constant for the four models. Note the bolometric luminosity of the starburst peaks 
at about 10-20Myrs. 



with the starburst galaxies. RRC proposed that the overlap is 
due to mixing of the starburst and Disc components. Rowan- 
Robinson (1992) also showed that some of the variation in the 
colours of the disc galaxies may be due to the intensity of the 
radiation field. 

If we consider first the predicted colours of the sequence of 
GMCs (dotted line), we see that they span the entire range of 
observed galaxy colours. This suggests that a weighted sum of 
emission from such a family of GMCs (which is what equation 
(12) essentially is) may explain the observed galaxy colours. To 
test this we have computed the colours of a galaxy that experi- 
enced an exponentially decaying burst with r = 20Myrs (solid 
line). The predicted colours in the age range 0-72Myrs nicely 
match the spread in the colours of starburst galaxies. Further- 
more, they predict a correlation of the 100/60 ratio, and an an- 
ticorrelation of the 25/12 one, with age. The predicted mean 
age of starburst galaxies (~ 50Myrs) agrees remarkably well 
with other estimates from Br7 equivalent widths and CO in- 
dices (Goldader et al 1997) as well as ISO spectroscopy (Genzel 
et al 1998). The predicted colours are not very sensitive on the 
assumed value of r. The colours for the constant star forma- 
tion case, a scenario more appropriate for disc galaxies, also lie 
on the same track but are packed towards the M82 end. This 
may have significant implications for the origin of the far-IR 
luminosity of disc galaxies as we discuss further in section 8. 
The conclusion from this analysis is therefore that the age of 
the burst can account for some of the variation in IRAS galaxy 
colours attributed by RRC to mixing with a cirrus component. 

There are two more features of the colours of the sequence 
of GMCs that are worth mentioning. First the discontinuity at 
t — 20Mvrs is due to the switch from the HII to supernova 



the right of the colour-colour diagrams will depend on the de- 
tails of the supernova phase which we treat very crudely in 
this paper. If the supernova ejecta expel the dust further away 
than we have assumed, then (in the absence of any other heat- 
ing source) we would expect the dust to cool down further and 
the 60/25 and 100/60 ratios to increase. The discontinuity in 
the colours may also be an artifact of our simplified treatment 
of the GMC evolution. In particular, stellar winds from young 
stars may impart as much kinetic energy in the ISM as the 
type II supernovae. Secondly, the apparent 'retrograde' motion 
of the tracks, best seen in the 25/12 versus 60/25 colour-colour 
diagram, is due to the radiative transfer effect, highlighted in 
section 3, which leads to the temporary heating up of the dust 
shell despite the fact that the radiation field is diminishing. 
The reality of these features can only be assessed by more de- 
tailed simulations of the evolution of GMCs. 

In Figure 3 we give the SEDs of the starburst discussed 
above (r = 20Myrs) for some representative ages. No correc- 
tion for leaking starlight is applied to the models. A number 
of features are worth highlighting: (1) there is a general ten- 
dency for the peak of the SED to shift to longer wavelengths 
with age. (2) The PAH features get stronger with age. (3) the 
9.7/im silicate feature gets shallower with age. (4) while the 
optical/UV light of the youngest bursts is almost completely 
obscured, the oldest starbursts are predicted to emit signifi- 
cantly in the optical/UV. This illustrates the extent to which 
old diffuse (and therefore optically thin) GMCs dominate the 
emission of old bursts. 
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6.1 Observational constraints on M82 

M82 being the nearest example of a starburst (D=3.25Mpc 
Tammann & Sandage 1968) has received a lot of attention 
both observationally and theoretically. The luminosity of the 
starburst (3 x 10 10 Lq Telesco & Harper 1980) is sufficient to 
outshine any pre-existing stellar population in the nucleus of 
this dwarf galaxy. The IR luminosity arises mostly in the cen- 
tral 200x400pc region (13" x 26") which shows a bilobal dis- 
tribution and is aligned with a central stellar bar. The bar 
is mainly made up of an old stellar population (Telesco et al 
1991) and it is thought that it has probably formed about 10 s 
years ago during the interaction of M81 with M82. Radio in- 
terferometry suggests a supernova frequency of 1 every 10-20 
years in the central region which implies a high star formation 
rate (Muxlow et al 1994). 

Rieke et al (1980), in their classic study of M82, drew 
attention to the fact that the starburst emission is spatially 
correlated at 10-, 20/xm and radio wavelengths. The emission 
region constitutes an elongated structure aligned with the ma- 
jor axis of the galaxy. The broad correspondence of the mid-IR 
emission with the radio was also confirmed by Telesco & Gezari 
(1992). The map of Telesco & Gezari shows the two prominent 
peaks seen in earlier observations but their higher spatial res- 
olution also allowed them to resolve a number of clumps with 
characteristic radii of about 20pc. This provides further confir- 
mation of the basic assumption of our model that the starburst 
is made of an ensemble of massive GMCs. 

Rieke et al also pointed out that the starburst shows a 
completely different appearance in the K band with the 'nu- 
cleus' of the galaxy being the prominent feature. There is no 
10[im feature associated with the nucleus. In the K band there 
is also extended emission along a stellar bar about a kpc in 
extent. Rieke et al also found that 2/im spectra of the nucleus 
show the absorption bands characteristic of giants and super- 
giants. Their interpretation of these findings is that the nucleus 
represents the first generation of stars formed about 5 x 10 7 yrs 
ago in an (exponentially decaying) starburst which has since 
propagated outwards. 

Shen & Lo (1995) mapped the CO emission from M82 
with 2.5" resolution. Their maps revealed unresolved (< 30pc) 
structure in the CO emission. They suggest that most of the 
CO emitting gas (and therefore the starburst) is located in 
molecular spilar arms at 125 and 390pc from the nucleus. They 
also resolved the previously known double-peaked structure in 
several peaks. 

The task of utilizing the multiwavelength data for M82 
is complicated by the fact that the good spatial and spectral 
resolution data at the shorter wavelengths A < 60/xm can not 
be exactly matched by the far-IR and submm data. A scan by 
Telesco & Harper (1980) at 58/im (where the spectral energy 
distribution peaks) shows that the intrinsic source size is about 
30" very similar to that at 10/xm. Hughes et al (1994) decon- 
volved the source size of M82 at various wavelengths. They find 
that between 10-40/wn the source size is fairly constant (25x8 
arcsec 2 ) but between 100pm and 1.3mm it is slightly larger 
(37x10 arcsec 2 ). 

While an IRAS LRS spectrum (with an aperture large 
enough to include the whole galaxy) is available for M82, be- 
cause of the importance of the mid and near-IR wavelength 
range in constraining the model, we use instead the 2-13/im 
spectrophotometry by Willner et al (1977) with a 30" aper- 
ture. The 10/j.m spectrum looks similar to that of Roche et al 
(1991) and Jones & Rodriguez- Espinosa (1984) in a number of 
nositions with 3" which lends further support to the idea that 



16 




log (V/im) 

Figure 4. Fit to the spectrum of M82 with the starburst model. As 
discussed in the text the model assumes that M82 has experienced 
two exponentially decaying bursts over the last 26.5Myrs. The first 
one (26. 5Myrs, r = 2Myrs) is shown by the dotted line and the 
second one (16.5Mj/rs, r = 6Myrs) by a dashed line. The total 
emission is given by the solid line. Data from Willner et al 1977, 
Johnson et al (1966), RRE and references therein 

forming regions. The IRAS LRS spectrum also shows a similar 
spectrum although the 10/xm feature appears a little shallower. 
The flux in the IRAS LRS spectrum is also about a factor of 2 
higher than the spectrum of Willner but also about 50% higher 
than the broad band 12pm point. The 10pm growth curves of 
Rieke et al (1980) also suggest that the total mid-IR flux may 
be higher than that in the 30" aperture by about 50%. The K 
band photometry of Rieke in a 30" aperture matches that of 
Willner et al in a similar aperture rather well. The K flux in the 
30" aperture is about 30% of the total (Rieke et al, Ichikawa 
et al 1995). 

In conclusion, the data plotted in Figure ^ are probably 
a good description of the spectrum of the starburst in M82 
but it should be borne in mind that the far-IR and submm 
data may include contributions from a region more extended 
than the central starburst. Having experimented with a range 
of combinations of age and decay time r we could not obtain 
a fit to the multi- wavelength data with a single burst scenario. 
The fundamental problem with a single burst (in the context of 
our model) is that M82 shows characteristics for both a young 
burst (warm SED and fairly deep 10pm absorption) and an 
older burst (strong 1-5/xm continuum). 

Our model fit, shown in Figure ^, assumes that M82 
experienced two bursts of star formation. The first one oc- 
cured 26.5Myrs ago and had a very steep exponential decay 
(r = 2Myrs). This burst is now responsible for most of the 
K band light and about half of the far-IR and submillimeter 
emission. The second burst occured 16.5Myrs ago and decayed 
more slowly (r = GMyrs). The GMCs in both bursts are as- 
sumed to leak 20% of their starlight after lOMyrs which subse- 
quently suffers a visual extinction of 1.5 magnitudes. The two 
bursts are predicted to contribute roughly equal amounts of 
UV flux. The spectrophotometry of Willner et al shows an ex- 
cess over the model longwards of the 3.3pm feature (possibly 
associated with a quasi-continuum similar to that underlying 
the longer wavelength features) which we cannot match with 
our present grain model. Note that while the 10pm spectrum 
of the younger burst provides an excellent fit to the shape of 
the observed spectrum, the addition of the older burst with 
its flatter spectrum dilutes the absorption feature. Our model 
could be in better agreement with the 10«m feature and the 
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Figure 5. Fit to the spectrum of NGC6090 with an exponentially 
decaying starburst model (t = 64Myrs,T = 50Myrs). Data from 
Acosta-Pulido, J. A., et al, (1996), and Gordon et al (1997) 



shells associated with the older burst have expanded beyond 
the central 30 arcsec. In fact there is evidence from the map of 
Hughes et al (1994) that some submm emission is associated 
with the outflows along the minor axis of the galaxy. 



7 A MODEL FOR NGC6090 

Acosta-Pulido et al (1996) presented ISO spectrophotome- 
try and extensive multi-wavelength photometry for NGC6090. 
This galaxy is about 10 times more luminous than M82 and 
its optical image shows a disturbed morphology and signs of a 
recent interaction. UV to K band data for this galaxy are also 
available (Gordon et al 1997). The IRAS colours of NGC6090 
show evidence for cold dust. RRC attributed that to a signifi- 
cant contribution from cirrus in this galaxy. It is therefore not 
surprising that the emission from an M82 type starburst has to 
be supplemented by colder dust to fit the data for A > 120/^m 
(Acosta-Pulido et al). 

As it is clear from Figure 2, the IRAS colours of NGC6090 
predict an 'old' burst scenario. In Figure ^ we test this by com- 
paring such a model to the multiwavelength data. We find that 
a good fit can be obtained with t — 64Myrs, r = 50Myrs. As 
in the M82 model we have assumed that after lOMyrs 20% of 
the starlight from each GMC escapes (and subsequently suffers 
1.5mag of visual extinction), but we find that adequate fits to 
the data can be obtained even if we have no leakage. 

The remaining discrepancy between the model and the far- 
IR observations suggests that some emission from even colder 
dust in the starburst or general interstellar medium of this 
galaxy (cirrus) may be contributing at the longer wavelengths. 
Colder dust in the starburst could for example arise if the su- 
pernova explosions are more effective (as we discuss in section 
5) in destroying the GMCs and expelling the dust away from 
the heating sources. Evidence of very cold dust (T ~ 10K) has 
been found in some normal galaxies by recent ISO observations 
(e.g. Kriigel et al 1998). The same authors, however, find no 
evidence for cold dust in the starburst galaxies in their sample. 
It will be interesting to assess how common is this far-IR ex- 
cess in starburst galaxies as it will provide strong constraints 
to theoretical models. 



8 SUMMARY AND DISCUSSION 

We have described an evolutionary scheme for massive GMCs 
centrally illuminated by a stellar cluster. The defining charac- 
teristic of this scheme is that the HII regions formed by the 
ionizing stellar radiation compress the neutral gas and dust 
in the GMC to a narrow shell. This naturally explains why 
the maximum temperature of the large grains in HII regions 
appears to be lower than the sublimation temperature of the 
graphite and silicate grains. It also explains why the near to 
mid-IR spectrum of HII regions and starburst galaxies is dom- 
inated by the PAH emission. These trends are consistent with 
the findings of Helou et al (1998) from ISO mid-IR spectropho- 
tometry. The evolution of the GMC after about 10 7 yrs is more 
uncertain but we show that even a modest rate of expansion 
at the interstellar velocity dispersion is sufficient to lead to 
the formation of a diffuse cold shell by about 10 8 yrs. The in- 
frared spectrum of these shells is very similar to infrared cirrus 
clouds. We associate these shells with the HI superstructures 
observed in our galaxy and other galaxies (Heiles 1979) and 
usually associated with large star clusters. 

The sequence of GMC spectral energy distributions we 
have computed cover the range of observed galaxy SEDs as 
observed by IRAS. This leads us to believe that the observed 
spectra of galaxies could be modelled in terms of this evolu- 
tionary scheme and constrain their star formation history. Ex- 
ponentially decaying bursts are shown to account satisfactorily 
for the colours of starburst galaxies in the sample of RRC. The 
age of the burst is shown to be an important contributing fac- 
tor to the spread in the galaxy colours and can mimick (to some 
extent) the effect of mixing with a cirrus component. Our ap- 
plication of these models to M82 and NGC6090 has shown that 
good fits to the UV to mm data on these galaxies can be ob- 
tained with the models. In M82 we find evidence for two bursts 
separated by lOMyrs. Our model is similar in this respect to 
that of Rieke et al (1993) who find that two Gaussian bursts 
(separated by 8-25Myrs) can account for a range of observa- 
tional constraints. Theoretical support for this periodic burst 
scenario is also provided by the models of Kriigel & Tutukov 
(1993). NGC6090 is best fitted by an older burst (64Myrs). 
There is evidence for a colder dust component in this galaxy. 
It is not clear whether that is due to cirrus or colder dust in 
the starburst. 

The colours predicted for a continuous star formation sce- 
nario, one more appropriate for disc galaxies, falls short of 
explaining the colours of disc galaxies. Although a more thor- 
ough exploration of the parameter space is needed before any 
definite conclusions can be drawn, this result is not entirely sur- 
prising for the following reason. In starburst galaxies, almost 
by definition, the emission from the recently formed stars out- 
shines the old stellar population so the effect of the latter on 
our models is probably not very significant. In disc galaxies, 
however, where the star formation rate is lower, the old stel- 
lar population is more important and needs to be taken into 
account when considering the energy balance and emission of 
GMCs, especially the older ones. The effect of this would be, 
neglecting any change in the spectral shape, to boost the lumi- 
nosity of the old GMCs and therefore their contribution to the 
overall galaxy emission. This provides an indirect argument 
in favour of the idea that old stars make a significant contri- 
bution to the far-IR luminosity of disc galaxies (Waltebros & 
Greenwalt 1996). 

The wealth of data on M82 allows us to apply further 
checks on the validity of this model. O'Connell et al (1995) 
imaged the central few hundred parsecs of M82 with the HST 
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Figure 6. Predicted evolution of the luminosity (normalized to its 
peak) of a starburst with r = 20Myrs at different infrared wave- 
lengths. 



luminous star clusters. The clusters do not seem to be associ- 
ated with X-rays, infrared or radio compact sources but instead 
with the less obscured regions in the galaxy such as the 'nu- 
cleus' and the periphery of the dust lane. The size of these 
clusters is estimated to be about 3pc, consistent with being 
the remnants of the GMCs in our model for the starburst. The 
fact that M82 is viewed almost edge-on complicates somewhat 
the interpretation of these observations but the (at least) par- 
tial segregation of the optical and infrared emitting regions 
suggests that they arise from stellar populations at different 
evolutionary stages, as in our model. This also suggests that 
the assumption of spherical symmetry is reasonable at least for 
the young [t < lOMyrs) GMCs. 

Satyapal et al (1997) studied the stellar clusters in the 
central 500pc of M82 with near-IR spectroscopy and high res- 
olution imaging. Their findings support a picture in which the 
typical age of the stellar clusters is 10 7 yrs. They also find a cor- 
relation between the age of the clusters and distance from the 
centre suggesting that the starburst is propagating outwards 
at a speed of about 50km/s. This picture is also supported 
by spectro-imaging of M82 at 3.3/im which shows evidence for 
dissociation of PAH molecules in the 'nucleus' (Normand et 
al 1995). McLeod et al (1993) also found that the 3.3/xra fea- 
ture is about twice as strong at a position offset by about 8" 
from the nucleus. The centre of M82 has been known for some 
time (Axon & Taylor 1978) to show evidence for a biconical 
outflow along its minor axis. This is now widely believed to 
be due to a galactic wind, a general characteristic of starburst 
galaxies, powered by stellar winds and supernovae. Numerical 
hydrodynamic simulations (Suchkov et al 1994) suggest that 
such outflows can develop under quite general conditions at the 
center of starbursts but their morphology depends very much 
on the energy deposition rate and density of the medium. It is 
likely that the galactic wind was powered by the earlier burst 
and that the interaction of the wind with gas in the plane of 
the galaxy triggered or at least contributed to the outwardly 
propagating secondary burst. 

Telesco & Gezari (1992) noted that there is no correspon- 
dence between the peaks of the radio emission in M82 (pre- 
sumed to be young supernova remnants) and those at 12.4/ira. 
This suggests that the mid-IR peaks in M82 are associated 
with GMCs in which supernova activity has not yet started or 
in which the supernova rate is still very low. This observation 
is in excellent agreement with our M82 model in which the 
mid-IR emission is predominantly due to the vouneer burst. 



galaxy is predicted to change quite considerably with age. We 
illustrate this more clearly in Figure 6 where we plot the lumi- 
nosity at different IRAS wavelengths (normalized to its peak) 
as a function of age for the same model parameters used in Fig- 
ure 3 (r = 20Myrs). The luminosity rises sharply in the first 
lOMyrs to peak at an age which depends on the wavelength. 
It then declines sharply by up to a factor of 5 with age. We 
have omitted the 25/im curve as it is very similar to that at 
60/im. There are a number of interesting implications arising 
from this result for the detectability of starburst galaxies. It 
implies for example that mid-IR surveys would preferentially 
detect younger starbursts than far-IR surveys. It also implies 
that near the detection limit of such surveys up to half the 
galaxies that experienced a burst in the last 7 x 10 7 years may 
be missed because they are either too young or too old. This 
may have implications for estimates of the star formation rate 
over volumes of space as it introduces another form of bias. 

In conclusion, the evolutionary scheme we have put for- 
ward in this paper seems to be in very good agreement with an 
array of observational evidence on M82 and is consistent with 
the infrared properties of other starburst galaxies. It promises 
to be useful for the interpretation of the growing datasets on 
infrared galaxies and we plan such studies in future work. The 
scheme also lends itself for incorporation into simulations of 
galaxy formation and evolution and this was one of the prime 
motivations for its development. 

While useful constraints on the star formation history and 
stellar populations of galaxies can be obtained from this model, 
it will be of much interest to explore in the future the ef- 
fect of deviations from spherical symmetry especially at the 
later stages of GMC evolution. Hydrodynamical simulations 
of the evolution of star-forming molecular clouds which take 
into account the effects of ionization, stellar winds and multi- 
supernova explosions in dense environments should also take 
high priority. Only then will we be able to take full advantage 
of the data ISO, SCUBA, WIRE, SIRTF, NGST, VLT, FIRST, 
PLANCK, SOFIA, ALMA etc. are expected to yield over the 
next decade or so. 
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APPENDIX: CORRECTIONS FOR 
NON-SPHERICAL GEOMETRY 

Let the monochromatic luminosity of a spherical GMC of op- 
tical depth Ty (centrally illuminated by a stellar cluster) at 
frequency v be 4-^5^ whereas that of the stellar cluster (in 
the absence of any dust) be A-kSZ- By definition, J Sidv = 
J o °° Stdu. If the optical depth of the GMC is t v (< t v ) for 
a fraction / of the sky centred on the star cluster, then the 
luminosity of such a non-spherical GMC can be approximated 
by, 

L„=4n[(l-f)SZ+fe v St+f J ° St] 

Sidv 

J o 

where C e ,v is the dust extinction cross-section at frequency v. 

The first term assumes that the GMC emits isotropically 
with the same spectral energy distribution as the spherical 
GMC. The factor 1 — /, however, accounts for the fact that 
a fraction / of the sky is covered by a lower (or even zero) 
optical depth. The remaining two terms attempt to correct 
for this complication. The second term gives the stelar light 
that is transmitted through the optically thin holes. The third 
term assumes that the light absorbed by dust in the optically 
thin holes is reprocessed to the infrared with the same spectral 
energy distribution as that of the spherical GMC. 

If a large number of non-spherical, but otherwise identi- 
cal, GMCs viewed at different orientations are present in a 
starburst, then in equations 12 and 13 one can use their av- 
erage emission S v = While this approximation conserves 
energy and is probably close to the best one can do for cor- 
recting for non-spherical geometry, without doing the actual 
calculation, its usefulness is bound to be limited as some of 
the assumptions (e.g. isotropic emission from a non-spherical 
optically thick system) are known to be false (e.g. Efstathiou 
& Rowan- Robinson 1995). 



